Part I: Spatial Data

1 Types of spatial data

Spatial data are used across a wide range of fields to support decision-making, including environment, public health, ecology, agriculture, urban planning, economy, and society.

1.1 Areal data

Areal data usually arise when the number of events corresponding to some variable of interest are aggregated in areas.

library(sf)
Warning: package 'sf' was built under R version 4.4.3
Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE
library(mapview)
d <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
mapview(d, zcol = "SID74")

Example of areal data. Number of sudden infant deaths in counties of North Carolina, USA, in 1974.

library(spData)
Warning: package 'spData' was built under R version 4.4.2
library(ggplot2)
d <- st_read(system.file("shapes/columbus.gpkg", package = "spData"), quiet = TRUE)
st_crs(d) <- NA
ggplot(d) + geom_sf(aes(fill = INC))

Example of areal data. Household income in $1000 USD in neighborhoods in Columbus, Ohio, in 1980.

library(terra)
Warning: package 'terra' was built under R version 4.4.3
terra 1.8.60
d <- rast(system.file("ex/elev.tif", package = "terra"))
plot(d)

Example of areal data. Elevation at raster grid cells covering Luxembourg.

1.2 Geostatistical data

For example, we can use air pollution measurements at a number of monitoring stations to predict air pollution at other locations taking into account spatial autocorrelation and other factors that are known to predict the outcome of interest.

library(sp)
library(sf)
library(mapview)
data(meuse)
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 28992)
mapview(meuse, zcol = "lead", map.types = "CartoDB.Voyager")

Example of geostatistical data. Topsoil lead concentrations at locations sampled in a flood plain of the river Meuse, The Netherlands.

library(spData)
mapview(properties, zcol = "prpsqm")

Example of geostatistical data. Price per square meter of a set of apartments in Athens, Greece, in 2017.

library(malariaAtlas)
Warning: package 'malariaAtlas' was built under R version 4.4.3
Registered S3 method overwritten by 'malariaAtlas':
  method           from   
  autoplot.default ggplot2
d <- getPR(country = "Zimbabwe", species = "BOTH")
Loading ISO 19139 XML schemas...
Loading ISO 19115-3 XML schemas...
Loading ISO 19139 codelists...
Please Note: Because you did not provide a version, by default the version being used is 202406 (This is the most recent version of PR data. To see other version options use function listPRPointVersions)
Warning: The `sourcedata` argument of `isAvailable_pr()` is deprecated as of
malariaAtlas 1.6.0.
ℹ The argument 'sourcedata' has been deprecated. It will be removed in the next
  version. It has no meaning.
ℹ The deprecated feature was likely used in the malariaAtlas package.
  Please report the issue at
  <https://github.com/malaria-atlas-project/malariaAtlas/issues>.
Please Note: Because you did not provide a version, by default the version being used is 202406 (This is the most recent version of PR data. To see other version options use function listPRPointVersions)
Confirming availability of PR data for: Zimbabwe...
PR points are available for Zimbabwe.
Attempting to download PR point data for Zimbabwe ...
Data downloaded for Zimbabwe.
NOTE: All available data for this query was downloaded for both species,
 but there are no PR points for P. vivax in this region in the MAP database. 
To check endemicity patterns or to contribute data, visit malariaatlas.org OR email us at map@bdi.ox.ac.uk.
ggplot2::autoplot(d)
Warning: The `format` argument of `getShp()` is deprecated as of malariaAtlas 1.6.0.
ℹ The argument 'format' has been deprecated. It will be removed in the next
  version. Admin boundaries will be correctly plotted using autoplot without
  the argument.
ℹ The deprecated feature was likely used in the malariaAtlas package.
  Please report the issue at
  <https://github.com/malaria-atlas-project/malariaAtlas/issues>.
Please Note: Because you did not provide a version, by default the version being used is 202403 (This is the most recent version of admin unit shape data. To see other version options use function listShpVersions)

Example of geostatistical data. Malaria prevalence at specific locations in Zimbabwe.

1.3 Point patterns

Point patterns arise when the variable to be analyzed corresponds to the location of events.

library(spatstat)
Warning: package 'spatstat' was built under R version 4.4.2
Loading required package: spatstat.data
Warning: package 'spatstat.data' was built under R version 4.4.2
Loading required package: spatstat.univar
Warning: package 'spatstat.univar' was built under R version 4.4.2
spatstat.univar 3.1-1
Loading required package: spatstat.geom
Warning: package 'spatstat.geom' was built under R version 4.4.2
spatstat.geom 3.3-5

Attaching package: 'spatstat.geom'
The following objects are masked from 'package:terra':

    area, delaunay, is.empty, rescale, rotate, shift, where.max,
    where.min
Loading required package: spatstat.random
Warning: package 'spatstat.random' was built under R version 4.4.2
spatstat.random 3.3-2
Loading required package: spatstat.explore
Warning: package 'spatstat.explore' was built under R version 4.4.2
Loading required package: nlme
spatstat.explore 3.3-4
Loading required package: spatstat.model
Warning: package 'spatstat.model' was built under R version 4.4.2
Loading required package: rpart
spatstat.model 3.3-4
Loading required package: spatstat.linnet
Warning: package 'spatstat.linnet' was built under R version 4.4.2
spatstat.linnet 3.2-5

spatstat 3.3-1 
For an introduction to spatstat, type 'beginner' 
plot(clmfires, use.marks = FALSE, pch = ".")

Data clmfires is a marked point pattern containing information of each fire, depicts the location of the fires without the mark.

library(spatstat)
plot(hamster)

This figure also shows the positions of cell nuclei in a histological section of a tissue from a lymphoma in the kidney of a hamster from spatstat. The nuclei are classified as either “pyknotic” (corresponding to dying cells) or “dividing” (corresponding to cells arrested in the act of dividing).

library(sparr)
Warning: package 'sparr' was built under R version 4.4.3


Welcome to
   _____ ___  ____  ____  ____         
  / ___// _ \/ _  \/ __ \/ __ \        
  \__ \/ ___/ __  /  ___/  ___/        
 ___/ / /  / / / / /\ \/ /\ \          
/____/_/  /_/ /_/_/  \__/  \_\   v2.3-16

- type news(package="sparr") for an overview
- type help("sparr") for documentation
- type citation("sparr") for how to cite
data(pbc)
plot(unmark(pbc[which(pbc$marks == "case"), ]), main = "cases")
axis(1)
axis(2)
title(xlab = "Easting", ylab = "Northing")

plot(unmark(pbc[which(pbc$marks == "control"), ]), pch = 3, main = "controls")
axis(1)
axis(2)
title(xlab = "Easting", ylab = "Northing")

Examples of point patterns. Top: Locations of fires in Castilla-La Mancha, Spain, between 1998 and 2007. Bottom: Locations and types of cells in a tissue.

1.4 Spatio-temporal data

Spatio-temporal data arise when information is both spatially and temporally referenced. Thus, we can consider spatial data as temporal aggregations or temporal snapshots of a spatio-temporal process.

1.5 Spatial functional data

Spatial functional data arise when the three types of spatial data (areal, geostatistical, and point patterns) are combined with random functions.

1.6 Mobility data

Besides the three classical types of spatial data (i.e., areal, geostatistical, and point patterns), we can also consider other spatial data such as flows containing the number of individuals or other elements moving between locations.

library("epiflows")
Warning: package 'epiflows' was built under R version 4.4.3
epiflows is loaded with the following global variables in `global_vars()`:
coordinates, pop_size, duration_stay, first_date, last_date, num_cases
data("Brazil_epiflows")

loc <- merge(x = YF_locations, y = YF_coordinates,
             by.x = "location_code", by.y = "id", sort = FALSE)
ef <- make_epiflows(flows = YF_flows, locations = loc,
                   coordinates = c("lon", "lat"),
                   pop_size = "location_population",
                   duration_stay = "length_of_stay",
                   num_cases = "num_cases_time_window",
                   first_date = "first_date_cases",
                   last_date = "last_date_cases")
ef

/// Epidemiological Flows //

  // class: epiflows, epicontacts
  // 15 locations; 100 flows; directed
  // optional variables: coordinates, pop_size, duration_stay, num_cases, first_date, last_date 

  // locations
Warning: `as_data_frame()` was deprecated in tibble 2.0.0.
ℹ Please use `as_tibble()` (with slightly different semantics) to convert to a
  tibble, or `as.data.frame()` to convert to a data frame.
ℹ The deprecated feature was likely used in the epiflows package.
  Please report the issue at <https://github.com/reconhub/epiflows/issues>.
# A tibble: 15 × 8
   id                 location_population num_cases_time_window first_date_cases
   <chr>                            <dbl>                 <dbl> <fct>           
 1 Espirito Santo                 3973697                  2600 2017-01-04      
 2 Minas Gerais                  20997560                  4870 2016-12-19      
 3 Rio de Janeiro                16635996                   170 2017-02-19      
 4 Sao Paulo                     44749699                   200 2016-12-17      
 5 Southeast Brazil              86356952                  7840 2016-12-17      
 6 Argentina                           NA                    NA <NA>            
 7 Chile                               NA                    NA <NA>            
 8 Germany                             NA                    NA <NA>            
 9 Italy                               NA                    NA <NA>            
10 Paraguay                            NA                    NA <NA>            
11 Portugal                            NA                    NA <NA>            
12 Spain                               NA                    NA <NA>            
13 United Kingdom                      NA                    NA <NA>            
14 United States of …                  NA                    NA <NA>            
15 Uruguay                             NA                    NA <NA>            
# ℹ 4 more variables: last_date_cases <fct>, length_of_stay <dbl>, lon <dbl>,
#   lat <dbl>

  // flows

# A tibble: 100 × 3
   from             to         n
   <chr>            <chr>  <dbl>
 1 Espirito Santo   Italy  2828.
 2 Minas Gerais     Italy 15714.
 3 Rio de Janeiro   Italy  8164.
 4 Sao Paulo        Italy 34039.
 5 Southeast Brazil Italy 76282.
 6 Espirito Santo   Spain  3270.
 7 Minas Gerais     Spain 18176.
 8 Rio de Janeiro   Spain  9443.
 9 Sao Paulo        Spain 39371.
10 Southeast Brazil Spain 88231.
# ℹ 90 more rows

We can use this data to create an epiflows object called ef that allows us to use the prediction and visualization functions.

2 Spatial data in R

2.1 Vector data

The sf package allows us to work with vector data which is used to represent points, lines, and polygons.

2.1.1 Shapefile

Vector data are commonly stored in the shapefile format.
A shapefile is not a single file, but a collection of related files that work together.

The three mandatory components are:

  • .shp – contains the geometry data,

  • .shx – a positional index that allows quick access to the geometry,

  • .dbf – stores attribute data for each feature.

    Additional optional files may include:

  • .prj – specifies the map projection,

  • .sbn and .sbx – spatial index files,

  • .shp.xml – contains spatial metadata in XML format.

library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
map<-st_read(pathshp,quiet= TRUE)
class(map)
[1] "sf"         "data.frame"
head(map)
Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
  NWBIR74 BIR79 SID79 NWBIR79                       geometry
1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
2      10   542     3      12 MULTIPOLYGON (((-81.23989 3...
3     208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
4     123   830     2     145 MULTIPOLYGON (((-76.00897 3...
5    1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...
6     954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...
plot(map[1]) #plotfirstattribute

Map of the first attribute of the ’sf‘ object representing the counties of North Carolina, USA.

2.2 Raster data

Raster data (also referred to as grid data) is a spatial data structure that divides the region of study into rectangles of the same size called cells or pixels, and that can store one or more values for each of these cells.

2.2.1 GeoTIFF

Raster data often come in GeoTIFF format which has extension .tif.

library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- terra::rast(pathraster)
r
class       : SpatRaster 
size        : 90, 95, 1  (nrow, ncol, nlyr)
resolution  : 0.008333333, 0.008333333  (x, y)
extent      : 5.741667, 6.533333, 49.44167, 50.19167  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326) 
source      : elev.tif 
name        : elevation 
min value   :       141 
max value   :       547 
plot(r)

Map of raster data representing the elevation of Luxembourg.

2.3 Coordinate Reference System

2.3.1 Geographic CRS

In a geographic CRS, latitude and longitude values are used to identify locations on the Earth’s three-dimensional ellipsoid surface.

Latitude values range from –90° at the South Pole to +90° at the North Pole. Longitude values measure the angle east or west of the Prime Meridian.

Latitude and longitude coordinates may be expressed in degrees, minutes, and seconds, or in decimal degrees.

Projected CRS

Projected CRSs use Cartesian coordinates to reference a location on a two-dimensional representation of the Earth.

World maps with Mercator (left) and Robinson (right) projections.

A common projection is the Universal Transverse Mercator (UTM) projection. This projection is conformal, meaning it preserves angles and therefore maintains shapes over small areas.

2.3.3 EPSG codes

Most common CRSs can be specified using their EPSG (European Petroleum Survey Group) codes or their Proj4 strings. Common spatial projections can be found at https://spatialreference.org/ref/. Details of a given projection can be inspected using the st_crs() function of the sf package. For example, EPSG code 4326 refers to the WGS84 longitude/latitude projection.

st_crs("EPSG:4326")$Name
[1] "WGS 84"
 st_crs("EPSG:4326")$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
st_crs("EPSG:4326")$epsg
[1] 4326

2.3.4 Transforming CRS with sf and terra

library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
map <- st_read(pathshp, quiet = TRUE)
# Get CRS
# st_crs(map)
# Transform CRS
map2 <- st_transform(map, crs = "EPSG:4326")
# Get CRS
# st_crs(map2)
library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
# Get CRS
# crs(r)
# Transform CRS
r2 <- terra::project(r, "EPSG:2169")
# Get CRS
# crs(r2)

Alternatively, if we want the transformed data to align exactly with other raster data we are using, we can project it using an existing raster that has the desired geometry.

# x is existing raster
# r is raster we project
# r2 <- terra::project(r, x)

2.4 Old spatial packages

Before the sf package was developed, the sp package was used to represent and work with vector spatial data. The sp, rgdal (Bivand et al., 2023), rgeos (Bivand and Rundel, 2022), and maptools (Bivand and Lewin-Koh, 2022) packages are no longer maintained and have been retired.

3 The sf package for spatial vector data

3.1 The sf package

library(sf)
pathshp <- system.file("shape/nc.shp", package = "sf")
nc <- st_read(pathshp, quiet = TRUE)
class(nc)
[1] "sf"         "data.frame"

A sf object contains the following objects of class sf, sfc and sfg: • sf (simple feature): each row of the data.frame is a single simple feature consisting of attributes and geometry. • sfc (simple feature geometry list-column): the geometry column of the data.frame is a list-column of class sfc with the geometry of each simple feature. • sfg (simple feature geometry): each of the rows of the sfc list-column corresponds to the simple feature geometry (sfg) of a single simple feature.

head(nc)
Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
  NWBIR74 BIR79 SID79 NWBIR79                       geometry
1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...
2      10   542     3      12 MULTIPOLYGON (((-81.23989 3...
3     208  3616     6     260 MULTIPOLYGON (((-80.45634 3...
4     123   830     2     145 MULTIPOLYGON (((-76.00897 3...
5    1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...
6     954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...
plot(nc)
Warning: plotting the first 10 out of 14 attributes; use max.plot = 14 to plot
all

object representing the counties of North Carolina, USA, and associated information.

nc[1,] # first row
Simple feature collection with 1 feature and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
1 0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
  BIR79 SID79 NWBIR79                       geometry
1  1364     0      19 MULTIPOLYGON (((-81.47276 3...
nc[nc$NAME == "Ashe", ] # row with NAME "Ashe"
Simple feature collection with 1 feature and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
Geodetic CRS:  NAD27
   AREA PERIMETER CNTY_ CNTY_ID NAME  FIPS FIPSNO CRESS_ID BIR74 SID74 NWBIR74
1 0.114     1.442  1825    1825 Ashe 37009  37009        5  1091     1      10
  BIR79 SID79 NWBIR79                       geometry
1  1364     0      19 MULTIPOLYGON (((-81.47276 3...
nc[1, "NWBIR74"] # first row, column with name NWBIR74
Simple feature collection with 1 feature and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.23436 xmax: -81.23989 ymax: 36.58965
Geodetic CRS:  NAD27
  NWBIR74                       geometry
1      10 MULTIPOLYGON (((-81.47276 3...
nc[1, "NWBIR74", drop = TRUE] # drop geometry
[1] 10
attr(,"class")
[1] "numeric"
# Geometries printed in abbreviated form
st_geometry(nc)
Geometry set for 100 features 
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
Geodetic CRS:  NAD27
First 5 geometries:
MULTIPOLYGON (((-81.47276 36.23436, -81.54084 3...
MULTIPOLYGON (((-81.23989 36.36536, -81.24069 3...
MULTIPOLYGON (((-80.45634 36.24256, -80.47639 3...
MULTIPOLYGON (((-76.00897 36.3196, -76.01735 36...
MULTIPOLYGON (((-77.21767 36.24098, -77.23461 3...
# View complete geometry by selecting one
st_geometry(nc)[[1]]
MULTIPOLYGON (((-81.47276 36.23436, -81.54084 36.27251, -81.56198 36.27359, -81.63306 36.34069, -81.74107 36.39178, -81.69828 36.47178, -81.7028 36.51934, -81.67 36.58965, -81.3453 36.57286, -81.34754 36.53791, -81.32478 36.51368, -81.31332 36.4807, -81.26624 36.43721, -81.26284 36.40504, -81.24069 36.37942, -81.23989 36.36536, -81.26424 36.35241, -81.32899 36.3635, -81.36137 36.35316, -81.36569 36.33905, -81.35413 36.29972, -81.36745 36.2787, -81.40639 36.28505, -81.41233 36.26729, -81.43104 36.26072, -81.45289 36.23959, -81.47276 36.23436)))

3.2 Creating a sf object

Simple feature geometries (sfg objects) can be, for example, of type POINT (single point), MULTIPOINT (set of points), or POLYGON (polygon), and can be created using the functions st_point(), st_multipoint(), and st_polygon(), respectively. Here, we create an sf object containing two single points, a set of points, and a polygon, each with one attribute.

library(sf)
library(ggplot2)

# Single points (point as a vector)
p1_sfg <- st_point(c(2, 2))
p2_sfg <- st_point(c(2.5, 3))

# Set of points (points as a matrix)
p <- rbind(
  c(6, 2),
  c(6.1, 2.6),
  c(6.8, 2.5),
  c(6.2, 1.5),
  c(6.8, 1.8)
)
mp_sfg <- st_multipoint(p)

# Polygon: sequence of points forming a closed, non-self-intersecting ring
# The first ring denotes the exterior ring;
# zero or more subsequent rings denote holes within it
p1 <- rbind(c(10, 0), c(11, 0), c(13, 2),
            c(12, 4), c(11, 4), c(10, 0))
p2 <- rbind(c(11, 1), c(11, 2), c(12, 2), c(11, 1))
pol_sfg <- st_polygon(list(p1, p2))

# Create sf object
p_sfc <- st_sfc(p1_sfg, p2_sfg, mp_sfg, pol_sfg)
df <- data.frame(v1 = c("A", "B", "C", "D"))
p_sf <- st_sf(df, geometry = p_sfc)

# Plot single points, multipoint, and polygon
ggplot(p_sf) +
  geom_sf(aes(col = v1), size = 3) +
  theme_bw()

sf object representing two single points, a set of points, and a polygon, with one attribute.

3.3 st_*() functions

Common functions to manipulate sf objects include:
st_read() – reads an sf object
st_write() – writes an sf object
st_crs() – gets or sets a new coordinate reference system (CRS)
st_transform() – transforms data to a new CRS
st_intersection() – intersects sf objects
st_union() – combines several sf objects into one
st_simplify() – simplifies an sf object
st_coordinates() – retrieves coordinates of an sf object
st_as_sf() – converts another object type to an sf object

We can read a sf object as follows:

library(sf)
library(ggplot2)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))

We can inspect the first rows of the sf object map with head().

head(map)
Simple feature collection with 6 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
Geodetic CRS:  NAD27
# A tibble: 6 × 15
   AREA PERIMETER CNTY_ CNTY_ID NAME   FIPS  FIPSNO CRESS_ID BIR74 SID74 NWBIR74
  <dbl>     <dbl> <dbl>   <dbl> <chr>  <chr>  <dbl>    <int> <dbl> <dbl>   <dbl>
1 0.114      1.44  1825    1825 Ashe   37009  37009        5  1091     1      10
2 0.061      1.23  1827    1827 Alleg… 37005  37005        3   487     0      10
3 0.143      1.63  1828    1828 Surry  37171  37171       86  3188     5     208
4 0.07       2.97  1831    1831 Curri… 37053  37053       27   508     1     123
5 0.153      2.21  1832    1832 North… 37131  37131       66  1421     9    1066
6 0.097      1.67  1833    1833 Hertf… 37091  37091       46  1452     7     954
# ℹ 4 more variables: BIR79 <dbl>, SID79 <dbl>, NWBIR79 <dbl>,
#   geometry <MULTIPOLYGON [°]>

We can delete some of the polygons by taking a subset of the rows of map.

# Delete polygon
map <- map[-which(map$FIPS %in% c("37125", "37051")), ]
ggplot(map) + geom_sf(aes(fill = SID79))

We can use st_union() with argument by_feature = FALSE to combine all geometries together.

# Combine geometries
ggplot(st_union(map, by_feature = FALSE) %>%
         st_sf()) + geom_sf()
although coordinates are longitude/latitude, st_union assumes that they are
planar

The boundaries of a map can be simplified with the st_simplify() function.

 # Simplify
map_utm <- st_transform(map, 32617)
map_simplified <- st_simplify(map_utm, dTolerance = 10000)
ggplot(map_simplified) +
  geom_sf(fill = "lightblue", color = "gray40") +
  theme_minimal() +
  labs(title = "Simplified North Carolina Map (10 km tolerance)")

3.4 Transforming point data to a sf object

The st_as_sf() function allows us to convert a foreign object into an sf object. For example, we can have a data frame containing the coordinates of several locations along with their attributes, which we can then convert into an sf object.

library(sf)
library(mapview)
d <- data.frame(
  place = c("London", "Paris", "Madrid", "Rome"),
  long = c(-0.118092, 2.349014, -3.703339, 12.496366),
  lat = c(51.509865, 48.864716, 40.416729, 41.902782),
  value = c(200, 300, 400, 500))
class(d)
[1] "data.frame"
dsf <- st_as_sf(d, coords = c("long", "lat"))
st_crs(dsf) <- 4326
class(dsf)
[1] "sf"         "data.frame"
mapview(dsf)

3.5 Counting the number of points within polygons

We can use the st_intersects() function from sf to count the number of points within the polygons of an sf object. The returned object is a list containing the feature IDs that intersect each polygon. We can then use the lengths() function to calculate the number of points inside each feature.

library(sf)
library(ggplot2)

# Map with divisions (sf object)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))

# Generate random points within the polygons (sf geometry list-column sfc)
points <- st_sample(map, size = 100, exact = FALSE)
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Convert points (sfc) into an sf object so it can be used in ggplot easily
points_sf <- st_sf(geometry = points)

# Plot map and points
 ggplot()+geom_sf(data= map)+geom_sf(data=points)

Points within polygons.

# Intersection: which points fall in which polygons
inter <- st_intersects(map, points_sf)
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Add point count to each polygon
map$count <- lengths(inter)

# Plot number of points within each polygon
 ggplot(map)+geom_sf(aes(fill = count))

Number of points within polygons.

3.6 Identifying polygons containing points

library(sf)
library(ggplot2)
library(dplyr)

Attaching package: 'dplyr'
The following object is masked from 'package:nlme':

    collapse
The following objects are masked from 'package:terra':

    intersect, union
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
# Map with divisions (sf object)
map <- read_sf(system.file("shape/nc.shp", package = "sf"))

# Generate random points within the polygons
points <- st_sample(map, size = 3) %>%
  st_as_sf()
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Find which polygon each point falls inside
# (note: first argument = points, second = polygons)
inter <- st_intersects(points, map)
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
although coordinates are longitude/latitude, st_intersects assumes that they
are planar
# Add column with the name of the area containing each point
points$areaname <- map[unlist(inter), "NAME", drop = TRUE]

# Show points with corresponding area name
points
Simple feature collection with 3 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -79.0806 ymin: 34.94656 xmax: -76.03757 ymax: 35.96621
Geodetic CRS:  NAD27
                           x areaname
1  POINT (-79.0806 34.94656)     Hoke
2 POINT (-76.03757 35.91925)  Tyrrell
3 POINT (-78.32428 35.96621) Franklin
# Map
ggplot(map) +
  geom_sf() +
  geom_sf(data = points) +
  geom_sf_label(
    data = map[unlist(inter), ],
    aes(label = NAME),
    nudge_y = 0.2
  )
Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data

3.7 Joining map and data

Sometimes, a map and its corresponding data are available separately, and we may want to create an sf object that combines the map with its related data for easier manipulation and plotting. This can be done by joining the map and data using the left_join() function from the dplyr package.

library(rnaturalearth)
Warning: package 'rnaturalearth' was built under R version 4.4.3
map <- ne_countries(returnclass = "sf")
library(wbstats)
Warning: package 'wbstats' was built under R version 4.4.3
indicators <- wb_search(pattern = "pollution")
d <- wb_data(indicator = "EN.ATM.PM25.MC.M3", start_date = 2016, end_date = 2016)

Then, we use the left_join() function from dplyr to merge the map and data, specifying the variable(s) to join by using the by argument. In this example, we use the ISO3 standard country codes instead of country names, as names may differ between the map and the data frame. Figure 3.7 shows the resulting map created using ggplot2.

library(dplyr)
library(ggplot2)
library(viridis)
Loading required package: viridisLite
map1 <- left_join(map, d, by = c("iso_a3" = "iso3c"))
ggplot(map1) + geom_sf(aes(fill = EN.ATM.PM25.MC.M3)) +
  scale_fill_viridis() + labs(fill = "PM2.5") + theme_bw()

PM2.5 values in each of the world countries in 2016.

Note that when using left_join(), the resulting object’s class matches that of the first argument.
Thus, left_join(sf_object, data.frame_object) produces an sf object, while left_join(data.frame_object, sf_object) produces a data frame.

map1 <- left_join(map, d, by = c("iso_a3" = "iso3c"))
class(map1)
[1] "sf"         "data.frame"
d1 <- left_join(d, map, by = c("iso3c" = "iso_a3"))
class(d1)
[1] "tbl_df"     "tbl"        "data.frame"

4 The terra package for raster and vector data

4.1 Raster Data

library(terra)
pathraster <- system.file("ex/elev.tif", package = "terra")
r <- rast(pathraster)
plot(r)

We can also use rast() to create a SpatRaster object by specifying the number of columns, number of rows, and the minimum and maximum x and y values.

r <- rast(ncol = 10, nrow = 10, xmin = -150, xmax = -80, ymin = 20, ymax = 60)
r
class       : SpatRaster 
size        : 10, 10, 1  (nrow, ncol, nlyr)
resolution  : 7, 4  (x, y)
extent      : -150, -80, 20, 60  (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84) 
nrow(r) # number of rows
[1] 10
ncol(r) # number of columns
[1] 10
dim(r) # dimension
[1] 10 10  1
ncell(r) # number of cells
[1] 100
values(r) <- 1:ncell(r)
r2 <- r * r
s <- c(r, r2)
plot(s[[2]]) # layer 2

Many generic functions can be used to operate with rasters, as shown in the following examples.

plot(min(s))

plot(r + r + 10)

plot(round(r))

plot(r == 1)

4.2 Vector data

pathshp <- system.file("ex/lux.shp", package = "terra")
v <- vect(pathshp)
# Longitude and latitude values
long <- c(-0.118092, 2.349014, -3.703339, 12.496366)
lat <- c(51.509865, 48.864716, 40.416729, 41.902782)
longlat <- cbind(long, lat)

# CRS
crspoints <- "+proj=longlat +datum=WGS84"

# Attributes for points
d <- data.frame(
  place = c("London", "Paris", "Madrid", "Rome"),
  value = c(200, 300, 400, 500)
)

# SpatVector object
pts <- vect(longlat, atts = d, crs = crspoints)
pts
 class       : SpatVector 
 geometry    : points 
 dimensions  : 4, 2  (geometries, attributes)
 extent      : -3.703339, 12.49637, 40.41673, 51.50986  (xmin, xmax, ymin, ymax)
 coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
 names       :  place value
 type        :  <chr> <num>
 values      : London   200
                Paris   300
               Madrid   400
# Plot
plot(pts)

4.3 Cropping, masking, and aggregating raster data

library(terra)
r <- geodata::worldclim_country(country = "Spain", var = "tavg", res = 10, path = tempdir())
plot(r)

We can average the temperature raster data over months with the mean() function.

r <- mean(r)
plot(r)

We also download the map of Spain with the rnaturalearth package (South, 2017), and delete the Canary Islands region.

# Map
library(ggplot2)
map <- rnaturalearth::ne_states("Spain", returnclass = "sf")
map <- map[-which(map$region == "Canary Is."), ] # delete region
ggplot(map) + geom_sf()

We obtain the spatial extent of the map with terra:ext(). Then, we can use crop() to remove the part of the raster that is outside the spatial extent.

# Cropping
sextent <- terra::ext(map)
r <- terra::crop(r, sextent)
plot(r)

We can use the mask() function to convert all values outside the map to NA.

# Masking
r <- terra::mask(r, vect(map))
plot(r)

The aggregate() function in terra can be used to group raster cells and create a new raster with lower resolution (i.e., larger cells).

# Aggregating
r <- terra::aggregate(r, fact = 20, fun = "mean", na.rm = TRUE)
plot(r)

4.4 Extracting raster values at points

library(terra)
# Raster (SpatRaster)
r <- rast(system.file("ex/elev.tif", package = "terra"))
# Polygons (SpatVector)
v <- vect(system.file("ex/lux.shp", package = "terra"))

We use the terra functions centroids() to get the centroids of the division polygons, and crds() to obtain their coordinates.

points <- crds(centroids(v))
plot(r)
plot(v, add = TRUE)
points(points)

Elevation raster, and division and centroids of polygons in Luxembourg.

We can obtain raster values at points using extract(), with the first argument as the SpatRaster object and the second argument as the data frame containing the points.

# data frame with the coordinates
points <- as.data.frame(points)
valuesatpoints <- extract(r, points)
cbind(points, valuesatpoints)
          x        y ID elevation
1  6.009082 50.07064  1       444
2  6.127425 49.86614  2       295
3  5.886502 49.80014  3       382
4  6.165081 49.92886  4       404
5  5.914545 49.93892  5       414
6  6.378449 49.78511  6       320
7  6.311601 49.54569  7       193
8  6.346395 49.68742  8       228
9  5.963503 49.64159  9       313
10 6.023816 49.52331 10       282
11 6.167624 49.61815 11       328
12 6.113598 49.75744 12       221

4.5 Extracting and averaging raster values within polygons

We can also use extract() to obtain raster values from SpatRaster objects within polygons of class SpatVector.

# Extracted raster cells within each polygon
head(extract(r, v, na.rm = TRUE))
  ID elevation
1  1       547
2  1       485
3  1       497
4  1       515
5  1       515
6  1       515
# Extracted raster cells and percentage of area covered within each polygon
head(extract(r, v, na.rm = TRUE, weights = TRUE))
  ID elevation     weight
1  1        NA 0.04545454
2  1        NA 0.10909091
3  1       529 0.24545454
4  1       542 0.46363635
5  1       547 0.68181816
6  1       535 0.11818181
# Average raster values by polygon
v$avg <- extract(r, v, mean, na.rm = TRUE)$elevation

# Area-weighted average raster values by polygon (weights = TRUE)
v$weightedavg <- extract(r, v, mean, na.rm = TRUE, weights = TRUE)$elevation
library(ggplot2)
library(tidyterra)
Warning: package 'tidyterra' was built under R version 4.4.3
Registered S3 method overwritten by 'tidyterra':
  method              from        
  autoplot.SpatRaster malariaAtlas

Attaching package: 'tidyterra'
The following object is masked from 'package:stats':

    filter
# Plot average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = avg)) + scale_fill_terrain_c()

# Plot area-weighted average raster values within polygons
ggplot(data = v) + geom_spatvector(aes(fill = weightedavg)) + scale_fill_terrain_c()

Average and area-weighted average of elevation values in each of the divisions of Luxembourg.

5 Making maps with R

Maps allow us to effectively communicate spatial information. Here, we demonstrate how to create both static and interactive maps using several mapping packages, including ggplot2 (Wickham et al., 2022a), leaflet (Cheng et al., 2022a), mapview (Appelhans et al., 2022), and tmap (Tennekes, 2022).

library(sf)
nameshp <- system.file("shape/nc.shp", package = "sf")
d <- st_read(nameshp, quiet = TRUE)
d$vble <- d$SID74
d$vble2 <- d$SID79

5.1 ggplot2

library(ggplot2)
library(viridis)
ggplot(d) + geom_sf(aes(fill = vble)) +
  scale_fill_viridis() + theme_bw()

Alternatively, we can specify a device driver (e.g., png, pdf), print the plot, and then close the device using dev.off().

#png("plot.png")
#ggplot(d) + geom_sf(aes(fill = vble)) + scale_fill_viridis() + theme_bw()
#dev.off()

5.2 leaflet

library(plotly)

Attaching package: 'plotly'
The following object is masked from 'package:ggplot2':

    last_plot
The following object is masked from 'package:stats':

    filter
The following object is masked from 'package:graphics':

    layout
g <- ggplot(d) + geom_sf(aes(fill = vble))
ggplotly(g)

The sf object passed to leaflet() must have a geographic coordinate reference system (CRS) using latitude and longitude (EPSG:4326). Here, we use the st_transform() function from sf to convert the data d, which has CRS EPSG:4267, to CRS EPSG:4326.

st_crs(d)$epsg
[1] 4267
d <- st_transform(d, 4326)
library(leaflet)
Warning: package 'leaflet' was built under R version 4.4.2
pal <- colorNumeric(palette = "YlOrRd", domain = d$vble)
l <- leaflet(d) %>% addTiles() %>%
  addPolygons(color = "white", fillColor = ~ pal(vble),
              fillOpacity = 0.8) %>%
  addLegend(pal = pal, values = ~vble, opacity = 0.8)
l

We can also use the addMiniMap() function to include an inset map.

l %>% addMiniMap()

For example, to save the leaflet map l as .png we can proceed as follows:

# webshot::install_phantomjs()
# Saves map.html
# library(htmlwidgets)
# saveWidget(widget = l, file = "map.html")
# Takes a screenshot of the map.html created and saves it as map.png
# library(webshot)
# webshot::install_phantomjs()
# webshot(url = "map.html", file = "map.png")

5.2 mapview

library(mapview)
mapview(d, zcol = "vble")

Maps created with mapview can also be customized to include elements such as legends and background maps.

library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
mapview(d, zcol = "vble", map.types = "CartoDB.DarkMatter",
        col.regions = pal, layer.name = "SDI")

An inset map can also be added by using the addMiniMap() function of leaflet.

map1 <- mapview(d, zcol = "vble")
leaflet::addMiniMap(map1@map)

5.4 Side-by-side plots with mapview

library(leaflet.extras2)
Warning: package 'leaflet.extras2' was built under R version 4.4.3
library(RColorBrewer)

pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))

# common legend
at <- seq(min(c(d$vble, d$vble2)), max(c(d$vble, d$vble2)), length.out = 8)
m1 <- mapview(d, zcol = "vble", map.types = "CartoDB.Positron",
              col.regions = pal, at = at)
m2 <- mapview(d, zcol = "vble2", map.types = "CartoDB.Positron",
              col.regions = pal, at = at)
m1 | m2

5.5 Synchronized maps with leaf sync

library(RColorBrewer)

pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))

# common legend
at <- seq(min(c(d$vble, d$vble2)), max(c(d$vble, d$vble2)), length.out = 8)

m1 <- mapview(d, zcol = "vble", map.types = "CartoDB.Positron", col.regions = pal, at = at)
m2 <- mapview(d, zcol = "vble2", map.types = "CartoDB.Positron", col.regions = pal, at = at)

m <- leafsync::sync(m1, m2)
m

5.6 tmap

The tmap package (Tennekes, 2022) allows us to create static and interactive maps composed of multiple shapes and layers, with various styles.

library(tmap)
Warning: package 'tmap' was built under R version 4.4.2
tmap_mode("plot")
ℹ tmap mode set to "plot".
tm_shape(d) + tm_polygons("vble")

5.7 Maps of point data

library(maps)
Warning: package 'maps' was built under R version 4.4.2

Attaching package: 'maps'
The following object is masked from 'package:viridis':

    unemp
d <- world.cities
# Select South Africa
d <- d[which(d$country.etc == "South Africa"), ]
# Transform data to sf object
d <- st_as_sf(d, coords = c("long", "lat"))
# Assign CRS
st_crs(d) <- 4326
d$vble <- d$pop
d$size <- sqrt(d$vble)/100

ggplot(d) + geom_sf(aes(col = vble, size = size)) +
  scale_color_viridis()

A leaflet map can be created using addCircles() specifying the radius and color of the points.

pal <- colorNumeric(palette = "viridis", domain = d$vble)

leaflet(d) %>%
  addTiles() %>%
  addCircles(
    lng = st_coordinates(d)[, 1],
    lat = st_coordinates(d)[, 2],
    radius = ~sqrt(vble) * 10,
    color = ~pal(vble),
    popup = ~name
  ) %>%
  addLegend(
    pal = pal,
    values = ~vble,
    position = "bottomright"
  )

To create a map with mapview, we set the size of the points with cex = “size”.

d$size <- sqrt(d$vble)
mapview(d, zcol = "vble", cex = "size")

Finally, we use tm_dots() to create a map with tmap.

tmap_mode("view")
ℹ tmap mode set to "view".
tm_shape(d) + tm_dots("vble", scale = sqrt(d$vble)/500, palette = "viridis")

── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_dots()`: migrate the argument(s) related to the scale of the
visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
tm_scale(<HERE>).[v3->v4] `tm_tm_dots()`: migrate the argument(s) related to the scale of the
visual variable `size` namely 'scale' (rename to 'values.scale') to size.scale
= tm_scale_continuous(<HERE>).
ℹ For small multiples, specify a 'tm_scale_' for each multiple, and put them in
  a list: 'size.scale = list(<scale1>, <scale2>, ...)'

5.8 Maps of raster data

library(terra)
filename <- system.file("ex/elev.tif", package = "terra")
r <- rast(filename)

# Transform data to sf object
d <- st_as_sf(as.data.frame(r, xy = TRUE), coords = c("x", "y"))
# Assign CRS
st_crs(d) <- 4326
# Plot
ggplot(d) + geom_sf() +
  geom_raster(data = as.data.frame(r, xy = TRUE),
              aes(x = x, y = y, fill = elevation))

To use the leaflet and mapview packages, we transform the data from class terra to RasterLayer with the raster::brick() function.

library(raster)

Attaching package: 'raster'
The following object is masked from 'package:plotly':

    select
The following object is masked from 'package:tidyterra':

    select
The following object is masked from 'package:dplyr':

    select
The following object is masked from 'package:nlme':

    getData
rb <- raster::brick(r)
pal <- colorNumeric("YlOrRd", values(r), na.color = "transparent")

leaflet() %>%
  addTiles() %>%
  addRasterImage(rb, colors = pal, opacity = 0.8) %>%
  addLegend(
    pal = pal,
    values = values(r),
    title = "Elevation"
  )

We create a map with mapview setting the title with the argument layer as follows:

mapview(rb, layer = "elevation")

We create a map with tmap by using tm_raster().

tm_shape(r) + tm_raster(title = "elevation", palette = "YlOrRd")
── tmap v3 code detected ───────────────────────────────────────────────────────
[v3->v4] `tm_tm_raster()`: migrate the argument(s) related to the scale of the
visual variable `col` namely 'palette' (rename to 'values') to col.scale =
tm_scale(<HERE>).
[v3->v4] `tm_raster()`: migrate the argument(s) related to the legend of the
visual variable `col` namely 'title' to 'col.legend = tm_legend(<HERE>)'
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "YlOrRd" is named
"brewer.yl_or_rd"

5.9 Mobility flows with flowmapblue

#devtools::install_github("FlowmapBlue/flowmapblue.R")
library(flowmapblue)

locations <- data.frame(
 id = c(1, 2, 3),
 name = c("New York", "London", "Rio de Janeiro"),
 lat = c(40.713543, 51.507425,-22.906241),
 lon = c(-74.011219,-0.127738,-43.180244)
)

The data frame flows has the number of people moving between origin and destination locations.

flows <- data.frame(
  origin = c(1, 2, 3, 2, 1, 3),
  dest = c(2, 1, 1, 3, 3, 2),
  count = c(42, 51, 50, 40, 22, 42)
)

Finally, we call the flowmapblue() function, passing the locations and flows data frames, along with the Mapbox access token, and specifying options such as clustering and animation. If mapboxAccessToken is set to NULL, the flows between locations will still be visualized, but without a background map.

#flowmapblue(locations, flows, mapboxAccessToken = NULL, clustering = TRUE, darkMode = TRUE, animation = FALSE)

6 R packages to download open spatial data

6.1 Administrative boundaries of countries

# devtools::install_github("ropensci/rnaturalearthhires")
library(rnaturalearth)
library(sf)
library(ggplot2)
library(viridis)
library(patchwork)

Attaching package: 'patchwork'
The following object is masked from 'package:raster':

    area
The following object is masked from 'package:spatstat.geom':

    area
The following object is masked from 'package:terra':

    area
map1 <- ne_countries(type = "countries", country = "Germany", scale = "medium", returnclass = "sf")
map2 <- rnaturalearth::ne_states("Germany", returnclass = "sf")
p1 <- ggplot(map1) + geom_sf()
p2 <- ggplot(map2) + geom_sf()
p1 + p2

6.2 Climatic data

library(geodata)
Warning: package 'geodata' was built under R version 4.4.2
d <- worldclim_country(country = "Jamaica", var = "tmin", path = tempdir())
terra::plot(mean(d), plg = list(title = "Min. temperature (C)"))

6.3 Precipitation

library("chirps")
Warning: package 'chirps' was built under R version 4.4.2
location <- data.frame(long = 100.523186, lat = 13.736717)
d <- get_chirps(location, dates = c("2020-01-01", "2022-12-31"), server = "ClimateSERV")
writing: substituting ENGCRS["Undefined Cartesian SRS with unknown unit"] for missing CRS
Fetching data from ClimateSERV 
Getting your request...
ggplot(d, aes(x = date, y = chirps)) + geom_line() +
  labs(y = "Precipitation (mm)")

6.4 Elevation

library(rnaturalearth)
library(elevatr)
Warning: package 'elevatr' was built under R version 4.4.3
elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
will be dropped in future versions, so please plan accordingly.
library(terra)

map <- ne_countries(type = "countries", country = "Switzerland", scale = "medium", returnclass = "sf")
d <- get_elev_raster(locations = map, z = 9, clip = "locations")
Mosaicing & Projecting
Clipping DEM to locations
Note: Elevation units are in meters.
terra::plot(rast(d), plg = list(title = "Elevation (m)"))

6.5 OpenStreetMap data

library(osmdata)
Warning: package 'osmdata' was built under R version 4.4.3
Data (c) OpenStreetMap contributors, ODbL 1.0. https://www.openstreetmap.org/copyright
head(available_features())
[1] "4wd_only"  "abandoned" "abutters"  "access"    "addr"      "addr:*"   
head(available_tags("amenity"))
# A tibble: 6 × 2
  Key     Value          
  <chr>   <chr>          
1 amenity animal_boarding
2 amenity animal_breeding
3 amenity animal_shelter 
4 amenity animal_training
5 amenity arts_centre    
6 amenity atm            
placebb <- getbb("Barcelona")
placebb
        min       max
x  2.052498  2.228356
y 41.317035 41.467914

For example, we can obtain the hospitals in Barcelona by specifying its bounding box placebb and using add_osm_feature() with key = "amenity" and value = "hospital" as follows.

hospitals <- placebb %>% opq() %>%
  add_osm_feature(key = "amenity", value = "hospital") %>%
  osmdata_sf()

Motorways can be retrieved using key = "highway" and value = "motorway".

motorways <- placebb %>% opq() %>%
  add_osm_feature(key = "highway", value = "motorway") %>%
  osmdata_sf()
library(leaflet)
leaflet() %>% addTiles() %>%
  addPolylines(data = motorways$osm_lines, color = "black") %>%
  addPolygons(data = hospitals$osm_polygons,
              label = hospitals$osm_polygons$name)

6.6 World Bank data

library(wbstats)
indicators <- wb_search(pattern = "poverty|unemployment")
# print(indicators)
d <- wb_data(indicator = "MO.INDEX.HDEV.XQ", start_date = 2011, end_date = 2011)
print(head(d))
# A tibble: 6 × 9
  iso2c iso3c country            date MO.INDEX.HDEV.XQ unit  obs_status footnote
  <chr> <chr> <chr>             <dbl>            <dbl> <chr> <chr>      <chr>   
1 AO    AGO   Angola             2011             47.7 <NA>  <NA>       <NA>    
2 BI    BDI   Burundi            2011             48.5 <NA>  <NA>       <NA>    
3 BJ    BEN   Benin              2011             53.2 <NA>  <NA>       <NA>    
4 BF    BFA   Burkina Faso       2011             45.8 <NA>  <NA>       <NA>    
5 BW    BWA   Botswana           2011             80.3 <NA>  <NA>       <NA>    
6 CF    CAF   Central African …  2011             32.9 <NA>  <NA>       <NA>    
# ℹ 1 more variable: last_updated <date>
library(rnaturalearth)
library(mapview)
map <- ne_countries(continent = "Africa", returnclass = "sf")
map <- dplyr::left_join(map, d, by = c("iso_a3" = "iso3c"))
mapview(map, zcol = "MO.INDEX.HDEV.XQ")

6.7 Species occurrence

library('spocc')
Warning: package 'spocc' was built under R version 4.4.3
df <- occ(query = "Bradypus variegatus", from = "gbif",
          date = c("2000-01-01", "2019-12-31"),
          gbifopts = list(country = "CR"),
          has_coords = TRUE, limit = 1000)
d <- occ2df(df)
library(sf)
d <- st_as_sf(d, coords = c("longitude", "latitude"))
st_crs(d) <- 4326
mapview(d)